function [yobsCounter,statesCounter]=...
    kfilterNanRegCounterfactual(parvec,funcmod,Y,solveopt,addsol,...
    etaOrig,aZero)
% =========================================================================

%% 1. Solve model 
[G,R,C,eu,SDX,Z]=feval(funcmod,parvec,solveopt,addsol);
if ~isequal(eu,[1;1]);return;end 

%% 2. Initial Dimensions 

%% 3.1. Dimensions and storage 
[NObs,NZ]=size(Y); 
NS=size(G,1); 
NX=size(SDX,1); 

[yobsCounter,statesCounter]=feval(@kfilterRegSplitSimulation,aZero,etaOrig');

   
%% Subroutine kfilterRegSplitSimulation allows to simulate the model  
% Inputs 
% sInitial: [NS 1] Initial State 
% innovMat: [Nx Nobx] matrix of innovations 
% By being a sub-routine it has access to all variables defined above 
% Be-careful not to use an index (ii,jj) used above or to repeat variable
% names
    function [ySim,sSim]=kfilterRegSplitSimulation(sInitial,innovMat)
        if ~isequal(size(innovMat),[NX NObs])
            error('Input innovMat must be [NX NObs]')
        end
        sSim=zeros(NS,NObs);
        ySim=zeros(NZ,NObs);
        sSim(:,1)=G*sInitial+R*innovMat(:,1);
        ySim(:,1)=Z*(sSim(:,1)+C);
        for kk=2:NObs;
            sSim(:,kk)=G*sSim(:,kk-1)+R*innovMat(:,kk);
            ySim(:,kk)=Z*(sSim(:,kk)+C);
        end
        ySim=ySim';
        sSim=sSim';
    end

%% End of File
end